Data integration - QC

suppressMessages(library(Matrix))
suppressMessages(library(SingleCellExperiment))
suppressMessages(library(biomaRt))
suppressMessages(library(scDblFinder))
suppressMessages(library(Seurat))
suppressMessages(library(ggplot2))
suppressMessages(library(scran))
suppressMessages(library(scater))
suppressMessages(library(ggpubr))

Upload sparse matrices:

matrix_dir1_mct1 = "/mnt/nmorais-nfs/marta/pA_karine/r_session/e12h-preliminary-analysis/cellranger/mct1/outs/filtered_feature_bc_matrix/"
barcode.path_mct1 <- paste0(matrix_dir1_mct1, "barcodes.tsv.gz")
features.path_mct1 <- paste0(matrix_dir1_mct1, "features.tsv.gz")
matrix.path_mct1 <- paste0(matrix_dir1_mct1, "matrix.mtx.gz")
mat_mct1 <- readMM(file = matrix.path_mct1)
feature.names_mct1 = read.delim(features.path_mct1, 
                           header = FALSE,
                           stringsAsFactors = FALSE)
barcode.names_mct1 = read.delim(barcode.path_mct1, 
                           header = FALSE,
                           stringsAsFactors = FALSE)

colnames(mat_mct1) = paste0(barcode.names_mct1$V1, "_mct1")
rownames(mat_mct1) = feature.names_mct1$V2
feature.ids_mct1 = feature.names_mct1$V1
mat_mct1[1:5,1:5]

## 5 x 5 sparse Matrix of class "dgTMatrix"
##         AAACCCAAGTGGAAGA-1_mct1 AAACCCACAGTGTGGA-1_mct1 AAACCCACATGTGCTA-1_mct1
## Xkr4                          .                       .                       .
## Gm1992                        .                       .                       .
## Gm19938                       .                       .                       .
## Gm37381                       .                       .                       .
## Rp1                           .                       .                       .
##         AAACCCAGTTGGTACT-1_mct1 AAACGAAAGATCACTC-1_mct1
## Xkr4                          .                       .
## Gm1992                        .                       .
## Gm19938                       .                       .
## Gm37381                       .                       .
## Rp1                           .                       .

matrix_dir3_mct3 = "/mnt/nmorais-nfs/marta/pA_karine/r_session/e12h-preliminary-analysis/cellranger/mct3/outs/filtered_feature_bc_matrix/"
barcode.path_mct3 <- paste0(matrix_dir3_mct3, "barcodes.tsv.gz")
features.path_mct3 <- paste0(matrix_dir3_mct3, "features.tsv.gz")
matrix.path_mct3 <- paste0(matrix_dir3_mct3, "matrix.mtx.gz")
mat_mct3 <- readMM(file = matrix.path_mct3)
feature.names_mct3 = read.delim(features.path_mct3, 
                           header = FALSE,
                           stringsAsFactors = FALSE)
barcode.names_mct3 = read.delim(barcode.path_mct3, 
                           header = FALSE,
                           stringsAsFactors = FALSE)

colnames(mat_mct3) = paste0(barcode.names_mct3$V1, "_mct3")
rownames(mat_mct3) = feature.names_mct3$V2
feature.ids_mct3 = feature.names_mct3$V1
mat_mct3[1:5,1:5]

## 5 x 5 sparse Matrix of class "dgTMatrix"
##         AAACGAACATAATCCG-1_mct3 AAACGAATCGTGAGAG-1_mct3 AAACGCTAGTGCTCAT-1_mct3
## Xkr4                          .                       .                       .
## Gm1992                        .                       .                       .
## Gm19938                       .                       .                       .
## Gm37381                       .                       .                       .
## Rp1                           .                       .                       .
##         AAACGCTAGTTCTCTT-1_mct3 AAACGCTCATTACGGT-1_mct3
## Xkr4                          .                       .
## Gm1992                        .                       .
## Gm19938                       .                       .
## Gm37381                       .                       .
## Rp1                           .                       .

matrix_dir1_nt1 = "/mnt/nmorais-nfs/marta/pA_karine/r_session/e12h-preliminary-analysis/cellranger/nt1/outs/filtered_feature_bc_matrix/"
barcode.path_nt1 <- paste0(matrix_dir1_nt1, "barcodes.tsv.gz")
features.path_nt1 <- paste0(matrix_dir1_nt1, "features.tsv.gz")
matrix.path_nt1 <- paste0(matrix_dir1_nt1, "matrix.mtx.gz")
mat_nt1 <- readMM(file = matrix.path_nt1)
feature.names_nt1 = read.delim(features.path_nt1, 
                           header = FALSE,
                           stringsAsFactors = FALSE)
barcode.names_nt1 = read.delim(barcode.path_nt1, 
                           header = FALSE,
                           stringsAsFactors = FALSE)

colnames(mat_nt1) = paste0(barcode.names_nt1$V1, "_nt1")
rownames(mat_nt1) = feature.names_nt1$V2
feature.ids_nt1 = feature.names_nt1$V1
mat_nt1[1:5,1:5]

## 5 x 5 sparse Matrix of class "dgTMatrix"
##         AAACCCAAGAAACTGT-1_nt1 AAACCCAAGAGACAAG-1_nt1 AAACCCAAGTAAGACT-1_nt1
## Xkr4                         2                      .                      .
## Gm1992                       .                      .                      .
## Gm19938                      .                      .                      .
## Gm37381                      .                      .                      .
## Rp1                          .                      .                      .
##         AAACCCAGTAATGTGA-1_nt1 AAACCCAGTGCTATTG-1_nt1
## Xkr4                         .                      .
## Gm1992                       .                      .
## Gm19938                      .                      .
## Gm37381                      .                      .
## Rp1                          .                      .

matrix_dir_s1 = "/mnt/nmorais-nfs/Backup_Unit/Datasets/KarineSerre/matrices/filtered_feature_bc_matrix_1/"
barcode.path_s1 <- paste0(matrix_dir_s1, "barcodes.tsv.gz")
features.path_s1 <- paste0(matrix_dir_s1, "features.tsv.gz")
matrix.path_s1 <- paste0(matrix_dir_s1, "matrix.mtx.gz")
mat_s1 <- readMM(file = matrix.path_s1)
feature.names_s1 = read.delim(features.path_s1, 
                           header = FALSE,
                           stringsAsFactors = FALSE)
barcode.names_s1 = read.delim(barcode.path_s1, 
                           header = FALSE,
                           stringsAsFactors = FALSE)

colnames(mat_s1) = paste0(barcode.names_s1$V1, "_s1")
rownames(mat_s1) = feature.names_s1$V2
mat_s1[1:5,1:5]

## 5 x 5 sparse Matrix of class "dgTMatrix"
##         AAACCCAGTAGCGAGT-1_s1 AAACGAAGTGACCTGC-1_s1 AAACGAAGTTGCGGCT-1_s1
## Xkr4                        .                     .                     .
## Gm1992                      .                     .                     .
## Gm37381                     .                     .                     .
## Rp1                         .                     .                     .
## Sox17                       .                     .                     .
##         AAACGAATCTCTGGTC-1_s1 AAAGAACAGTTGCGAG-1_s1
## Xkr4                        .                     .
## Gm1992                      .                     .
## Gm37381                     .                     .
## Rp1                         .                     .
## Sox17                       .                     .

matrix_dir_s2 = "/mnt/nmorais-nfs/Backup_Unit/Datasets/KarineSerre/matrices/filtered_feature_bc_matrix_2/"
barcode.path_s2 <- paste0(matrix_dir_s2, "barcodes.tsv.gz")
features.path_s2 <- paste0(matrix_dir_s2, "features.tsv.gz")
matrix.path_s2 <- paste0(matrix_dir_s2, "matrix.mtx.gz")
mat_s2 <- readMM(file = matrix.path_s2)
feature.names_s2 = read.delim(features.path_s2, 
                           header = FALSE,
                           stringsAsFactors = FALSE)
barcode.names_s2 = read.delim(barcode.path_s2, 
                           header = FALSE,
                           stringsAsFactors = FALSE)

colnames(mat_s2) = paste0(barcode.names_s2$V1, "_s2")
rownames(mat_s2) = feature.names_s2$V2
mat_s2[1:5,1:5]

## 5 x 5 sparse Matrix of class "dgTMatrix"
##         AAACCCAAGGGAACAA-1_s2 AAACGAACATCCGTGG-1_s2 AAACGAACATGGCACC-1_s2
## Xkr4                        .                     .                     .
## Gm1992                      .                     .                     .
## Gm37381                     .                     .                     .
## Rp1                         .                     .                     .
## Sox17                       .                     .                     .
##         AAACGAAGTAGCTGCC-1_s2 AAACGAAGTCAAAGTA-1_s2
## Xkr4                        .                     .
## Gm1992                      .                     .
## Gm37381                     .                     .
## Rp1                         .                     .
## Sox17                       .                     .

matrix_dir_s3 = "/mnt/nmorais-nfs/Backup_Unit/Datasets/KarineSerre/matrices/filtered_feature_bc_matrix_3/"
barcode.path_s3 <- paste0(matrix_dir_s3, "barcodes.tsv.gz")
features.path_s3 <- paste0(matrix_dir_s3, "features.tsv.gz")
matrix.path_s3 <- paste0(matrix_dir_s3, "matrix.mtx.gz")
mat_s3 <- readMM(file = matrix.path_s3)
feature.names_s3 = read.delim(features.path_s3, 
                           header = FALSE,
                           stringsAsFactors = FALSE)
barcode.names_s3 = read.delim(barcode.path_s3, 
                           header = FALSE,
                           stringsAsFactors = FALSE)

colnames(mat_s3) = paste0(barcode.names_s3$V1, "_s3")
rownames(mat_s3) = feature.names_s3$V2
mat_s3[1:5,1:5]

## 5 x 5 sparse Matrix of class "dgTMatrix"
##         AAACCCAAGACGGTCA-1_s3 AAACCCAAGCTAGCCC-1_s3 AAACCCACACAATTCG-1_s3
## Xkr4                        .                     .                     .
## Gm1992                      .                     .                     .
## Gm37381                     .                     .                     .
## Rp1                         .                     .                     .
## Sox17                       .                     .                     .
##         AAACCCACACGGGCTT-1_s3 AAACCCACAGCAAGAC-1_s3
## Xkr4                        .                     .
## Gm1992                      .                     .
## Gm37381                     .                     .
## Rp1                         .                     .
## Sox17                       .                     .

matrix_dir_s4 = "/mnt/nmorais-nfs/Backup_Unit/Datasets/KarineSerre/matrices/filtered_feature_bc_matrix_4/"
barcode.path_s4 <- paste0(matrix_dir_s4, "barcodes.tsv.gz")
features.path_s4 <- paste0(matrix_dir_s4, "features.tsv.gz")
matrix.path_s4 <- paste0(matrix_dir_s4, "matrix.mtx.gz")
mat_s4 <- readMM(file = matrix.path_s4)
feature.names_s4 = read.delim(features.path_s4, 
                           header = FALSE,
                           stringsAsFactors = FALSE)
barcode.names_s4 = read.delim(barcode.path_s4, 
                           header = FALSE,
                           stringsAsFactors = FALSE)

colnames(mat_s4) = paste0(barcode.names_s4$V1, "_s4")
rownames(mat_s4) = feature.names_s4$V2
mat_s4[1:5,1:5]

## 5 x 5 sparse Matrix of class "dgTMatrix"
##         AAACCCAAGGTTTGAA-1_s4 AAACCCACAATTGAAG-1_s4 AAACCCACAGCTTCGG-1_s4
## Xkr4                        .                     .                     .
## Gm1992                      .                     .                     .
## Gm37381                     .                     .                     .
## Rp1                         .                     .                     .
## Sox17                       .                     .                     .
##         AAACCCAGTGACATCT-1_s4 AAACCCATCGACGACC-1_s4
## Xkr4                        .                     .
## Gm1992                      .                     .
## Gm37381                     .                     .
## Rp1                         .                     .
## Sox17                       .                     .

matrix_dir_s5 = "/mnt/nmorais-nfs/Backup_Unit/Datasets/KarineSerre/matrices/filtered_feature_bc_matrix_5/"
barcode.path_s5 <- paste0(matrix_dir_s5, "barcodes.tsv.gz")
features.path_s5 <- paste0(matrix_dir_s5, "features.tsv.gz")
matrix.path_s5 <- paste0(matrix_dir_s5, "matrix.mtx.gz")
mat_s5 <- readMM(file = matrix.path_s5)
feature.names_s5 = read.delim(features.path_s5, 
                           header = FALSE,
                           stringsAsFactors = FALSE)
barcode.names_s5 = read.delim(barcode.path_s5, 
                           header = FALSE,
                           stringsAsFactors = FALSE)

colnames(mat_s5) = paste0(barcode.names_s5$V1, "_s5")
rownames(mat_s5) = feature.names_s5$V2
mat_s5[1:5,1:5]

## 5 x 5 sparse Matrix of class "dgTMatrix"
##         AAACCCAAGGAATGTT-1_s5 AAACCCAGTAACAGGC-1_s5 AAACCCAGTCGGCCTA-1_s5
## Xkr4                        .                     .                     .
## Gm1992                      .                     .                     .
## Gm37381                     .                     .                     .
## Rp1                         .                     .                     .
## Sox17                       .                     .                     .
##         AAACCCATCCGTAGGC-1_s5 AAACCCATCCTTTAGT-1_s5
## Xkr4                        .                     .
## Gm1992                      .                     .
## Gm37381                     .                     .
## Rp1                         .                     .
## Sox17                       .                     .

matrix_dir_s6 = "/mnt/nmorais-nfs/Backup_Unit/Datasets/KarineSerre/matrices/filtered_feature_bc_matrix_6/"
barcode.path_s6 <- paste0(matrix_dir_s6, "barcodes.tsv.gz")
features.path_s6 <- paste0(matrix_dir_s6, "features.tsv.gz")
matrix.path_s6 <- paste0(matrix_dir_s6, "matrix.mtx.gz")
mat_s6 <- readMM(file = matrix.path_s6)
feature.names_s6 = read.delim(features.path_s6, 
                           header = FALSE,
                           stringsAsFactors = FALSE)
barcode.names_s6 = read.delim(barcode.path_s6, 
                           header = FALSE,
                           stringsAsFactors = FALSE)

colnames(mat_s6) = paste0(barcode.names_s6$V1, "_s6")
rownames(mat_s6) = feature.names_s6$V2
mat_s6[1:5,1:5]

## 5 x 5 sparse Matrix of class "dgTMatrix"
##         AAACCCACAAGCGAAC-1_s6 AAACCCACATCCAATG-1_s6 AAACCCAGTTCAACGT-1_s6
## Xkr4                        .                     .                     .
## Gm1992                      .                     .                     .
## Gm37381                     .                     .                     .
## Rp1                         .                     .                     .
## Sox17                       .                     .                     .
##         AAACGAAAGATTGTGA-1_s6 AAACGAAGTCATATGC-1_s6
## Xkr4                        .                     .
## Gm1992                      .                     .
## Gm37381                     .                     .
## Rp1                         .                     .
## Sox17                       .                     .

matrix_dir_s7 = "/mnt/nmorais-nfs/Backup_Unit/Datasets/KarineSerre/matrices/filtered_feature_bc_matrix_7/"
barcode.path_s7 <- paste0(matrix_dir_s7, "barcodes.tsv.gz")
features.path_s7 <- paste0(matrix_dir_s7, "features.tsv.gz")
matrix.path_s7 <- paste0(matrix_dir_s7, "matrix.mtx.gz")
mat_s7 <- readMM(file = matrix.path_s7)
feature.names_s7 = read.delim(features.path_s7, 
                           header = FALSE,
                           stringsAsFactors = FALSE)
barcode.names_s7 = read.delim(barcode.path_s7, 
                           header = FALSE,
                           stringsAsFactors = FALSE)

colnames(mat_s7) = paste0(barcode.names_s7$V1, "_s7")
rownames(mat_s7) = feature.names_s7$V2
mat_s7[1:5,1:5]

## 5 x 5 sparse Matrix of class "dgTMatrix"
##         AAACCCACAAGGTTGG-1_s7 AAACCCACAGCTGTAT-1_s7 AAACGAAAGAAGCTCG-1_s7
## Xkr4                        .                     .                     .
## Gm1992                      .                     .                     .
## Gm37381                     .                     .                     .
## Rp1                         .                     .                     .
## Sox17                       .                     .                     .
##         AAACGAAAGGCTCTAT-1_s7 AAACGAACAAATGGCG-1_s7
## Xkr4                        .                     .
## Gm1992                      .                     .
## Gm37381                     .                     .
## Rp1                         .                     .
## Sox17                       .                     .

matrix_dir_s8 = "/mnt/nmorais-nfs/Backup_Unit/Datasets/KarineSerre/matrices/filtered_feature_bc_matrix_8/"
barcode.path_s8 <- paste0(matrix_dir_s8, "barcodes.tsv.gz")
features.path_s8 <- paste0(matrix_dir_s8, "features.tsv.gz")
matrix.path_s8 <- paste0(matrix_dir_s8, "matrix.mtx.gz")
mat_s8 <- readMM(file = matrix.path_s8)
feature.names_s8 = read.delim(features.path_s8, 
                           header = FALSE,
                           stringsAsFactors = FALSE)
barcode.names_s8 = read.delim(barcode.path_s8, 
                           header = FALSE,
                           stringsAsFactors = FALSE)

colnames(mat_s8) = paste0(barcode.names_s8$V1, "_s8")
rownames(mat_s8) = feature.names_s8$V2
mat_s8[1:5,1:5]

## 5 x 5 sparse Matrix of class "dgTMatrix"
##         AAACCCAAGTAAGACT-1_s8 AAACCCAAGTTCTACG-1_s8 AAACCCAGTAACTGCT-1_s8
## Xkr4                        .                     .                     .
## Gm1992                      .                     .                     .
## Gm37381                     .                     .                     .
## Rp1                         .                     .                     .
## Sox17                       .                     .                     .
##         AAACCCATCACGGGCT-1_s8 AAACCCATCTGGCCAG-1_s8
## Xkr4                        .                     .
## Gm1992                      .                     .
## Gm37381                     .                     .
## Rp1                         .                     .
## Sox17                       .                     .

matrix_dir_s9 = "/mnt/nmorais-nfs/Backup_Unit/Datasets/KarineSerre/matrices/filtered_feature_bc_matrix_9/"
barcode.path_s9 <- paste0(matrix_dir_s9, "barcodes.tsv.gz")
features.path_s9 <- paste0(matrix_dir_s9, "features.tsv.gz")
matrix.path_s9 <- paste0(matrix_dir_s9, "matrix.mtx.gz")
mat_s9 <- readMM(file = matrix.path_s9)
feature.names_s9 = read.delim(features.path_s9, 
                           header = FALSE,
                           stringsAsFactors = FALSE)
barcode.names_s9 = read.delim(barcode.path_s9, 
                           header = FALSE,
                           stringsAsFactors = FALSE)

colnames(mat_s9) = paste0(barcode.names_s9$V1, "_s9")
rownames(mat_s9) = feature.names_s9$V2
mat_s9[1:5,1:5]

## 5 x 5 sparse Matrix of class "dgTMatrix"
##         AAACCCACACCGAATT-1_s9 AAACCCAGTCTGTCAA-1_s9 AAACGAAAGACGCAGT-1_s9
## Xkr4                        .                     .                     .
## Gm1992                      .                     .                     .
## Gm37381                     .                     .                     .
## Rp1                         .                     .                     .
## Sox17                       .                     .                     .
##         AAACGAAAGACTCTAC-1_s9 AAACGAATCCCAGGCA-1_s9
## Xkr4                        .                     .
## Gm1992                      .                     .
## Gm37381                     .                     .
## Rp1                         .                     .
## Sox17                       .                     .

dim(mat_mct1)

## [1] 32285  4476

dim(mat_mct3)

## [1] 32285  4328

dim(mat_nt1)

## [1] 32285  3855

dim(mat_s1)

## [1] 31053  1923

dim(mat_s2)

## [1] 31053  1897

dim(mat_s3)

## [1] 31053  8500

dim(mat_s4)

## [1] 31053  6756

dim(mat_s5)

## [1] 31053  6191

dim(mat_s6)

## [1] 31053  1635

dim(mat_s7)

## [1] 31053  2120

dim(mat_s8)

## [1] 31053  4373

dim(mat_s9)

## [1] 31053  2880

new_exp <- cbind(as.matrix(mat_mct1),as.matrix(mat_mct3), as.matrix(mat_nt1))

## Warning in asMethod(object): sparse->dense coercion: allocating vector of size
## 1.1 GiB

## Warning in asMethod(object): sparse->dense coercion: allocating vector of size
## 1.0 GiB

old_exp <- cbind(as.matrix(mat_s1), as.matrix(mat_s2), as.matrix(mat_s3), as.matrix(mat_s4), as.matrix(mat_s5),
                 as.matrix(mat_s6), as.matrix(mat_s7), as.matrix(mat_s8), as.matrix(mat_s9))

## Warning in asMethod(object): sparse->dense coercion: allocating vector of size
## 2.0 GiB

## Warning in asMethod(object): sparse->dense coercion: allocating vector of size
## 1.6 GiB

## Warning in asMethod(object): sparse->dense coercion: allocating vector of size
## 1.4 GiB

## Warning in asMethod(object): sparse->dense coercion: allocating vector of size
## 1.0 GiB

final_mat <- merge(new_exp, old_exp, by = "row.names", all.x = TRUE, all.y = TRUE) #keep all genes

final_mat[1:5,1:5]

dim(final_mat)

## [1] 33502 48935

rownames(final_mat) <- final_mat$Row.names
final_mat <- final_mat[,-1]
final_mat[1:5,1:5]

final_mat[is.na(final_mat)] <- 0

filter_genes <- rowSums(final_mat)
keep_genes <- filter_genes > 0
table(keep_genes)

## keep_genes
## FALSE  TRUE 
##  7692 25810

final_mat <- final_mat[keep_genes,]
dim(final_mat)

## [1] 25810 48934

Create SingleCellExperiment object:

sce <- SingleCellExperiment(assays = list(counts = final_mat))

Define control features - mitochondrial genes

ensembl=useMart("ensembl")
ensembl = useDataset("mmusculus_gene_ensembl",mart=ensembl)
#ensembl_gene_id             
mt_genes <- getBM(attributes=c("ensembl_gene_id", "external_gene_name" ), filters="chromosome_name",values="MT",mart=ensembl)
mt_genes

rownames(sce) <- gsub("\\.", "-", rownames(final_mat))
grep("\\.", rownames(sce), value = TRUE)

## character(0)

mt_genes_sce <- rownames(sce)[rownames(sce) %in% mt_genes$external_gene_name]
mt_genes_sce

##  [1] "mt-Atp6" "mt-Atp8" "mt-Co1"  "mt-Co2"  "mt-Co3"  "mt-Cytb" "mt-Nd1" 
##  [8] "mt-Nd2"  "mt-Nd3"  "mt-Nd4"  "mt-Nd4l" "mt-Nd5"  "mt-Nd6"

sample <- c(rep("MCT1", dim(mat_mct1)[2]),
            rep("MCT3", dim(mat_mct3)[2]),
            rep("NT1",  dim(mat_nt1)[2]),
            rep("S1-NT", dim(mat_s1)[2]),
            rep("S2-MCT", dim(mat_s2)[2]),
            rep("S3-PT", dim(mat_s3)[2]),
            rep("S4-NT", dim(mat_s4)[2]),
            rep("S5-MCT", dim(mat_s5)[2]),
            rep("S6-PT", dim(mat_s6)[2]),
            rep("S7-PT", dim(mat_s7)[2]),
            rep("S8-NT", dim(mat_s8)[2]),
            rep("S9-MCT", dim(mat_s9)[2])
            )

table(sample)

## sample
##   MCT1   MCT3    NT1  S1-NT S2-MCT  S3-PT  S4-NT S5-MCT  S6-PT  S7-PT  S8-NT 
##   4476   4328   3855   1923   1897   8500   6756   6191   1635   2120   4373 
## S9-MCT 
##   2880

sorting <- c(rep("day4", dim(mat_mct1)[2]),
            rep("day4", dim(mat_mct3)[2]),
            rep("day4",  dim(mat_nt1)[2]),
            rep("day1", dim(mat_s1)[2]),
            rep("day1", dim(mat_s2)[2]),
            rep("day2", dim(mat_s3)[2]),
            rep("day2", dim(mat_s4)[2]),
            rep("day2", dim(mat_s5)[2]),
            rep("day3", dim(mat_s6)[2]),
            rep("day3", dim(mat_s7)[2]),
            rep("day3", dim(mat_s8)[2]),
            rep("day3", dim(mat_s9)[2])
            )

condition <- c(rep("MCT - 12h", dim(mat_mct1)[2]),
            rep("MCT - 12h", dim(mat_mct3)[2]),
            rep("NT",  dim(mat_nt1)[2]),
            rep("NT", dim(mat_s1)[2]),
            rep("MCT - 72h", dim(mat_s2)[2]),
            rep("PT", dim(mat_s3)[2]),
            rep("NT", dim(mat_s4)[2]),
            rep("MCT - 72h", dim(mat_s5)[2]),
            rep("PT", dim(mat_s6)[2]),
            rep("PT", dim(mat_s7)[2]),
            rep("NT", dim(mat_s8)[2]),
            rep("MCT - 72h", dim(mat_s9)[2])
            )

table(sample, condition)

##         condition
## sample   MCT - 12h MCT - 72h   NT   PT
##   MCT1        4476         0    0    0
##   MCT3        4328         0    0    0
##   NT1            0         0 3855    0
##   S1-NT          0         0 1923    0
##   S2-MCT         0      1897    0    0
##   S3-PT          0         0    0 8500
##   S4-NT          0         0 6756    0
##   S5-MCT         0      6191    0    0
##   S6-PT          0         0    0 1635
##   S7-PT          0         0    0 2120
##   S8-NT          0         0 4373    0
##   S9-MCT         0      2880    0    0

sce$sample <- sample 
sce$condition <- condition
sce$sorting <- sorting

lib_sizes <- colSums(counts(sce))
check_countMatrix <- counts(sce) > 0
total_features <- colSums(check_countMatrix)
  
mt_sizes <- colSums(counts(sce)[mt_genes_sce,])
mt_pct <- (mt_sizes/lib_sizes)*100
  
df <- data.frame(library_size = lib_sizes,
                   total_features = total_features,
                   percentage_mt_genes = mt_pct,
                 sample = sce$sample)

plot_histogram <- function(df_qc, x ,bins, title, xlab, min)
{ p<- gghistogram(data=df_qc, 
              x = x,
              y = "..count..",
              bins= bins,
              alpha = .2,
              color="black",
              fill="grey",
              title= title,
              ylab="#cells",
              xlab=xlab) + 
    theme_bw() +
    theme(text = element_text(size=20)) + 
    geom_vline(xintercept = min, linetype="dotted",
               color = "red", size=0.5) 
  
}

plot_scatterLibFeat <- function(df_qc, minx, miny)
{ 
  df_qc <- df_qc[order(df_qc$percentage_mt_genes, decreasing = FALSE),]

  p <- ggplot(data = df_qc, aes(log10_library_size, log10_total_features, color = percentage_mt_genes)) + 
    geom_point( size = 1, alpha = 1) + 
    scale_color_gradient2(midpoint = 12, low = "darkgreen", mid = "yellow", high="red") + 
    xlab("Library size") +
    ylab("Unique features") +
    labs(color = "") +
    theme(text = element_text(size=20), 
          panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
          panel.background = element_blank(), axis.line = element_line(colour = "black")) + 
     geom_vline(xintercept = minx, linetype="dotted",
               color = "red", size=0.5) + 
     geom_hline(yintercept = miny, linetype="dotted",
                color = "red", size=0.5) 
  
  
}

plot_scatterLibMt <- function(df_qc, minx)
{ 
  
  df_qc <- df_qc[order(df_qc$percentage_mt_genes, decreasing = FALSE),]
  
  
  p <- ggplot(data = df_qc, aes(log10_library_size, percentage_mt_genes, color = percentage_mt_genes)) + 
    geom_point( size = 1, alpha = 1) + 
    scale_color_gradient2(midpoint = 12, low = "darkgreen", mid = "yellow", high="red") + 
    xlab("Library size") +
    ylab("% MT-genes") +
    labs(color = "") +
    theme(text = element_text(size=20), 
          panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
          panel.background = element_blank(), axis.line = element_line(colour = "black")) + 
    geom_vline(xintercept = minx, linetype="dotted",
               color = "red", size=0.5) 
  
}

p <- plot_histogram(df_qc = df, 
               x = "library_size", 
               bins = 50, title = "library size",
               xlab = "total counts",
               min = 1000)

## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

df$log10_library_size <- log10(df$library_size)
df$log10_total_features <- log10(df$total_features)

g <- plot_histogram(df_qc = df, 
               x = "log10_library_size", 
               bins = 50, title = "library size",
               xlab = "log10(total counts)",
               min = 3)
p

g

g <- plot_histogram(df_qc = df[df$sample == "S1-NT",], 
               x = "log10_library_size", 
               bins = 50, title = "Library size: S1 - NT",
               xlab = "log10(total counts)",
               min = 3.5)
g

g <- plot_histogram(df_qc = df[df$sample == "S4-NT",], 
               x = "log10_library_size", 
               bins = 50, title = "Library size: S4 - NT",
               xlab = "log10(total counts)",
               min = 3.5)
g

g <- plot_histogram(df_qc = df[df$sample == "S8-NT",], 
               x = "log10_library_size", 
               bins = 50, title = "Library size: S8 - NT",
               xlab = "log10(total counts)",
               min = 3.5)
g

g <- plot_histogram(df_qc = df[df$sample == "S3-PT",], 
               x = "log10_library_size", 
               bins = 50, title = "Library size: S3 - PT",
               xlab = "log10(total counts)",
               min = 3.5)
g

g <- plot_histogram(df_qc = df[df$sample == "S6-PT",], 
               x = "log10_library_size", 
               bins = 50, title = "Library size: S6 - PT",
               xlab = "log10(total counts)",
               min = 3.5)
g

g <- plot_histogram(df_qc = df[df$sample == "S7-PT",], 
               x = "log10_library_size", 
               bins = 50, title = "Library size: S7 - PT",
               xlab = "log10(total counts)",
               min = 3.5)
g

g <- plot_histogram(df_qc = df[df$sample == "NT1",], 
               x = "log10_library_size", 
               bins = 50, title = "Library size: NT1",
               xlab = "log10(total counts)",
               min = 3.5)
g

g <- plot_histogram(df_qc = df[df$sample == "MCT1",], 
               x = "log10_library_size", 
               bins = 50, title = "Library size: MCT1",
               xlab = "log10(total counts)",
               min = 3.5)
g

g <- plot_histogram(df_qc = df[df$sample == "MCT3",], 
               x = "log10_library_size", 
               bins = 50, title = "Library size: MCT3",
               xlab = "log10(total counts)",
               min = 3.5)
g

g <- plot_histogram(df_qc = df[df$sample == "S2-MCT",], 
               x = "log10_library_size", 
               bins = 50, title = "Library size: S2 - MCT",
               xlab = "log10(total counts)",
               min = 3.5)
g

g <- plot_histogram(df_qc = df[df$sample == "S5-MCT",], 
               x = "log10_library_size", 
               bins = 50, title = "Library size: S5 - MCT",
               xlab = "log10(total counts)",
               min = 3.5)
g

g <- plot_histogram(df_qc = df[df$sample == "S9-MCT",], 
               x = "log10_library_size", 
               bins = 50, title = "Library size: S9 - MCT",
               xlab = "log10(total counts)",
               min = 3.5)
g

p <- plot_histogram(df_qc = df, 
               x = "total_features", 
               bins = 50, title = "Unique genes detected",
               xlab = "Total features",
               min = 1000)

df$log10_library_size <- log10(df$library_size)
df$log10_total_features <- log10(df$total_features)

g <- plot_histogram(df_qc = df, 
               x = "log10_total_features", 
               bins = 50, title = "Unique genes detected",
               xlab = "log10(total features)",
               min = 3)
p

g

g <- plot_histogram(df_qc = df[df$sample == "S1-NT",], 
               x = "log10_total_features", 
               bins = 50, title = "Unique genes detected: S1 - NT",
               xlab = "log10(total features)",
               min = 3)
g

g <- plot_histogram(df_qc = df[df$sample == "S4-NT",], 
               x = "log10_total_features", 
               bins = 50, title = "Unique genes detected: S4 - NT",
               xlab = "log10(total features)",
               min = 3)
g

g <- plot_histogram(df_qc = df[df$sample == "S8-NT",], 
               x = "log10_total_features", 
               bins = 50, title = "Unique genes detected: S8 - NT",
               xlab = "log10(total features)",
               min = 3)
g

g <- plot_histogram(df_qc = df[df$sample == "S3-PT",], 
               x = "log10_total_features", 
               bins = 50, title = "Unique genes detected: S3 - PT",
               xlab = "log10(total features)",
               min = 3)
g

g <- plot_histogram(df_qc = df[df$sample == "S6-PT",], 
               x = "log10_total_features", 
               bins = 50, title = "Unique genes detected: S6 - PT",
               xlab = "log10(total features)",
               min = 3)
g

g <- plot_histogram(df_qc = df[df$sample == "S7-PT",], 
               x = "log10_total_features", 
               bins = 50, title = "Unique genes detected: S7 - PT",
               xlab = "log10(total features)",
               min = 3)
g

g <- plot_histogram(df_qc = df[df$sample == "NT1",], 
               x = "log10_total_features", 
               bins = 50, title = "Unique genes detected: NT1",
               xlab = "log10(total features)",
               min = 3)
g

g <- plot_histogram(df_qc = df[df$sample == "MCT1",], 
               x = "log10_total_features", 
               bins = 50, title = "Unique genes detected: MCT1",
               xlab = "log10(total features)",
               min = 3)
g

g <- plot_histogram(df_qc = df[df$sample == "MCT3",], 
               x = "log10_total_features", 
               bins = 50, title = "Unique genes detected: MCT3",
               xlab = "log10(total features)",
               min = 3)
g

g <- plot_histogram(df_qc = df[df$sample == "S2-MCT",], 
               x = "log10_total_features", 
               bins = 50, title = "Unique genes detected: S2 - MCT",
               xlab = "log10(total features)",
               min = 3)
g

g <- plot_histogram(df_qc = df[df$sample == "S5-MCT",], 
               x = "log10_total_features", 
               bins = 50, title = "Unique genes detected: S5 - MCT",
               xlab = "log10(total features)",
               min = 3)
g

g <- plot_histogram(df_qc = df[df$sample == "S9-MCT",], 
               x = "log10_total_features", 
               bins = 50, title = "Unique genes detected: S9 - MCT",
               xlab = "log10(total features)",
               min = 3)
g

p <- plot_scatterLibFeat(df, 3.5, 3)
p

g <- plot_scatterLibMt(df, 3.5)
g + geom_hline(yintercept = 25)

p <- plot_scatterLibFeat(df[df$sample == "MCT1",], 3.5, 3) + ggtitle("MCT1")
p

p <- plot_scatterLibFeat(df[df$sample == "MCT3",], 3.5, 3)+ ggtitle("MCT3")
p

p <- plot_scatterLibFeat(df[df$sample == "NT1",], 3.5, 3)+ ggtitle("NT1")
p

p <- plot_scatterLibFeat(df[df$sample == "S1-NT",], 3.5, 3)+ ggtitle("S1-NT")
p

p <- plot_scatterLibFeat(df[df$sample == "S2-MCT",], 3.5, 3)+ ggtitle("S2-MCT")
p

p <- plot_scatterLibFeat(df[df$sample == "S3-PT",], 3.5, 3)+ ggtitle("S3-PT")
p

p <- plot_scatterLibFeat(df[df$sample == "S4-NT",], 3.5, 3)+ ggtitle("S4-NT")
p

p <- plot_scatterLibFeat(df[df$sample == "S5-MCT",], 3.5, 3)+ ggtitle("S5-MCT")
p

p <- plot_scatterLibFeat(df[df$sample == "S6-PT",], 3.5, 3)+ ggtitle("S6-PT")
p

p <- plot_scatterLibFeat(df[df$sample == "S7-PT",], 3.5, 3)+ ggtitle("S7-PT")
p

p <- plot_scatterLibFeat(df[df$sample == "S8-NT",], 3.5, 3)+ ggtitle("S8-NT")
p

p <- plot_scatterLibFeat(df[df$sample == "S9-MCT",], 3.5, 3)+ ggtitle("S9-MCT")
p

p <- plot_scatterLibMt(df[df$sample == "MCT1",], 3.5) + ggtitle("MCT1")
p

p <- plot_scatterLibMt(df[df$sample == "MCT3",], 3.5)+ ggtitle("MCT3")
p

p <- plot_scatterLibMt(df[df$sample == "NT1",], 3.5)+ ggtitle("NT1")
p

p <- plot_scatterLibMt(df[df$sample == "S1-NT",], 3.5)+ ggtitle("S1-NT")
p

p <- plot_scatterLibMt(df[df$sample == "S2-MCT",], 3.5)+ ggtitle("S2-MCT")
p

p <- plot_scatterLibMt(df[df$sample == "S3-PT",], 3.5)+ ggtitle("S3-PT")
p

p <- plot_scatterLibMt(df[df$sample == "S4-NT",], 3.5)+ ggtitle("S4-NT")
p

p <- plot_scatterLibMt(df[df$sample == "S5-MCT",], 3.5)+ ggtitle("S5-MCT")
p

p <- plot_scatterLibMt(df[df$sample == "S6-PT",], 3.5)+ ggtitle("S6-PT")
p

p <- plot_scatterLibMt(df[df$sample == "S7-PT",], 3.5)+ ggtitle("S7-PT")
p

p <- plot_scatterLibMt(df[df$sample == "S8-NT",], 3.5)+ ggtitle("S8-NT")
p

p <- plot_scatterLibMt(df[df$sample == "S9-MCT",], 3.5)+ ggtitle("S9-MCT")
p

median(df$percentage_mt_genes[df$sample == "MCT1"])

## [1] 1.214362

mean(df$percentage_mt_genes[df$sample == "MCT1"])

## [1] 3.760509

median(df$percentage_mt_genes[df$sample == "MCT3"])

## [1] 1.004837

mean(df$percentage_mt_genes[df$sample == "MCT3"])

## [1] 2.913248

median(df$percentage_mt_genes[df$sample == "NT1"])

## [1] 1.80533

mean(df$percentage_mt_genes[df$sample == "NT1"])

## [1] 6.208564

median(df$percentage_mt_genes[df$sample == "S1-NT"])

## [1] 2.775718

mean(df$percentage_mt_genes[df$sample == "S1-NT"])

## [1] 6.219134

median(df$percentage_mt_genes[df$sample == "S2-MCT"])

## [1] 4.166406

mean(df$percentage_mt_genes[df$sample == "S2-MCT"])

## [1] 14.95335

median(df$percentage_mt_genes[df$sample == "S3-PT"])

## [1] 2.961088

mean(df$percentage_mt_genes[df$sample == "S3-PT"])

## [1] 4.359633

median(df$percentage_mt_genes[df$sample == "S4-NT"])

## [1] 2.689016

mean(df$percentage_mt_genes[df$sample == "S4-NT"])

## [1] 5.426065

median(df$percentage_mt_genes[df$sample == "S5-MCT"])

## [1] 4.351457

mean(df$percentage_mt_genes[df$sample == "S5-MCT"])

## [1] 10.64891

median(df$percentage_mt_genes[df$sample == "S6-PT"])

## [1] 2.751902

mean(df$percentage_mt_genes[df$sample == "S6-PT"])

## [1] 4.709598

median(df$percentage_mt_genes[df$sample == "S7-PT"])

## [1] 2.769398

mean(df$percentage_mt_genes[df$sample == "S7-PT"])

## [1] 4.446721

median(df$percentage_mt_genes[df$sample == "S8-NT"])

## [1] 2.415459

mean(df$percentage_mt_genes[df$sample == "S8-NT"])

## [1] 4.517199

median(df$percentage_mt_genes[df$sample == "S9-MCT"])

## [1] 3.714046

mean(df$percentage_mt_genes[df$sample == "S9-MCT"])

## [1] 9.707103

df2 <- data.frame(sample = unique(df$sample),
                  condition = c("MCT - 12h", "MCT - 12h", "NT", "NT", "MCT - 72h", "PT", "NT", "MCT - 72h", "PT", "PT", "NT", "MCT - 72h"),
                  mean_mt = c(mean(df$percentage_mt_genes[df$sample == "MCT1"]),
                              mean(df$percentage_mt_genes[df$sample == "MCT3"]),
                              mean(df$percentage_mt_genes[df$sample == "NT1"]),
                              mean(df$percentage_mt_genes[df$sample == "S1-NT"]),
                              mean(df$percentage_mt_genes[df$sample == "S2-MCT"]),
                              mean(df$percentage_mt_genes[df$sample == "S3-PT"]),
                              mean(df$percentage_mt_genes[df$sample == "S4-NT"]),
                              mean(df$percentage_mt_genes[df$sample == "S5-MCT"]),
                              mean(df$percentage_mt_genes[df$sample == "S6-PT"]),
                              mean(df$percentage_mt_genes[df$sample == "S7-PT"]),
                              mean(df$percentage_mt_genes[df$sample == "S8-NT"]),
                              mean(df$percentage_mt_genes[df$sample == "S9-MCT"])
                  ),
                  median_mt = c(median(df$percentage_mt_genes[df$sample == "MCT1"]),
                              median(df$percentage_mt_genes[df$sample == "MCT3"]),
                              median(df$percentage_mt_genes[df$sample == "NT1"]),
                              median(df$percentage_mt_genes[df$sample == "S1-NT"]),
                              median(df$percentage_mt_genes[df$sample == "S2-MCT"]),
                              median(df$percentage_mt_genes[df$sample == "S3-PT"]),
                              median(df$percentage_mt_genes[df$sample == "S4-NT"]),
                              median(df$percentage_mt_genes[df$sample == "S5-MCT"]),
                              median(df$percentage_mt_genes[df$sample == "S6-PT"]),
                              median(df$percentage_mt_genes[df$sample == "S7-PT"]),
                              median(df$percentage_mt_genes[df$sample == "S8-NT"]),
                              median(df$percentage_mt_genes[df$sample == "S9-MCT"])))

df2$sample <- factor(df2$sample, levels = c("NT1","S1-NT", "S4-NT", "S8-NT",
                                            "S3-PT", "S6-PT", "S7-PT", 
                                            "MCT1","MCT3","S2-MCT", "S5-MCT", "S9-MCT"))
df2

ggplot(data=df2, aes(x=sample, y=mean_mt, fill=condition)) +
  geom_bar(stat="identity")+
  theme_minimal() +
  theme(axis.text.x = element_text(angle=45, hjust =1),
        text = element_text(size=20)) + labs(x = "Sample", y="Mean of % MT genes") +
scale_fill_manual(values= c("green3", "yellow3", "dodgerblue", "coral1"))

ggplot(data=df2, aes(x=sample, y=median_mt, fill=condition)) +
  geom_bar(stat="identity")+
  theme_minimal() +
  theme(axis.text.x = element_text(angle=45, hjust =1),
        text = element_text(size=20)) + labs(x = "Sample", y="Median of % MT genes") +
scale_fill_manual(values= c("green3", "yellow3", "dodgerblue", "coral1"))

keep_cells <- df$log10_library_size >= 3.5 & 
              df$log10_total_features >= 3 &
              df$percentage_mt_genes < 25
table(df$sample)

## 
##   MCT1   MCT3    NT1  S1-NT S2-MCT  S3-PT  S4-NT S5-MCT  S6-PT  S7-PT  S8-NT 
##   4476   4328   3855   1923   1897   8500   6756   6191   1635   2120   4373 
## S9-MCT 
##   2880

table(keep_cells, df$sample)

##           
## keep_cells MCT1 MCT3  NT1 S1-NT S2-MCT S3-PT S4-NT S5-MCT S6-PT S7-PT S8-NT
##      FALSE  668  462  565   296    432   863   964   1221   159   197   531
##      TRUE  3808 3866 3290  1627   1465  7637  5792   4970  1476  1923  3842
##           
## keep_cells S9-MCT
##      FALSE    404
##      TRUE    2476

total_cells <- table(df$sample)
names(total_cells) <- NULL
total_cells

##  [1] 4476 4328 3855 1923 1897 8500 6756 6191 1635 2120 4373 2880

discarded_cells <- table(keep_cells, df$sample)[1,]
discarded_cells

##   MCT1   MCT3    NT1  S1-NT S2-MCT  S3-PT  S4-NT S5-MCT  S6-PT  S7-PT  S8-NT 
##    668    462    565    296    432    863    964   1221    159    197    531 
## S9-MCT 
##    404

(discarded_cells/total_cells)*100

##  [1] 14.924039 10.674677 14.656291 15.392616 22.772799 10.152941 14.268798
##  [8] 19.722177  9.724771  9.292453 12.142694 14.027778

sce

## class: SingleCellExperiment 
## dim: 25810 48934 
## metadata(0):
## assays(1): counts
## rownames(25810): a A030001D20Rik ... Zzef1 Zzz3
## rowData names(0):
## colnames(48934): AAACCCAAGTGGAAGA-1_mct1 AAACCCACAGTGTGGA-1_mct1 ...
##   TTTGGTTGTAATTGGA-1_s9 TTTGTTGCAGTCTGGC-1_s9
## colData names(3): sample condition sorting
## reducedDimNames(0):
## mainExpName: NULL
## altExpNames(0):

sce <- sce[,keep_cells]
sce

## class: SingleCellExperiment 
## dim: 25810 42172 
## metadata(0):
## assays(1): counts
## rownames(25810): a A030001D20Rik ... Zzef1 Zzz3
## rowData names(0):
## colnames(42172): AAACCCAAGTGGAAGA-1_mct1 AAACCCACAGTGTGGA-1_mct1 ...
##   TTTGGTTGTAATTGGA-1_s9 TTTGTTGCAGTCTGGC-1_s9
## colData names(3): sample condition sorting
## reducedDimNames(0):
## mainExpName: NULL
## altExpNames(0):

table(rowSums(counts(sce)) > 0)

## 
## FALSE  TRUE 
##   127 25683

sce <- sce[rowSums(counts(sce)) > 0,]
sce

## class: SingleCellExperiment 
## dim: 25683 42172 
## metadata(0):
## assays(1): counts
## rownames(25683): a A030001D20Rik ... Zzef1 Zzz3
## rowData names(0):
## colnames(42172): AAACCCAAGTGGAAGA-1_mct1 AAACCCACAGTGTGGA-1_mct1 ...
##   TTTGGTTGTAATTGGA-1_s9 TTTGTTGCAGTCTGGC-1_s9
## colData names(3): sample condition sorting
## reducedDimNames(0):
## mainExpName: NULL
## altExpNames(0):

df <- df[keep_cells,]

saveRDS(sce, "/mnt/nmorais-nfs/marta/pA_karine/r_session/e12h-preliminary-analysis/integration/data/sce-total.rds")

saveRDS(df, "/mnt/nmorais-nfs/marta/pA_karine/r_session/e12h-preliminary-analysis/integration/data/qc-df.rds" )

sessionInfo()

## R version 4.1.2 (2021-11-01)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 18.04.6 LTS
## 
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.7.1
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.7.1
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
##  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
## [1] stats4    stats     graphics  grDevices utils     datasets  methods  
## [8] base     
## 
## other attached packages:
##  [1] ggpubr_0.5.0                scater_1.22.0              
##  [3] scran_1.22.1                scuttle_1.4.0              
##  [5] ggplot2_3.4.1               SeuratObject_4.1.3         
##  [7] Seurat_4.1.1                scDblFinder_1.8.0          
##  [9] biomaRt_2.50.3              SingleCellExperiment_1.16.0
## [11] SummarizedExperiment_1.24.0 Biobase_2.54.0             
## [13] GenomicRanges_1.46.1        GenomeInfoDb_1.30.1        
## [15] IRanges_2.28.0              S4Vectors_0.32.4           
## [17] BiocGenerics_0.40.0         MatrixGenerics_1.6.0       
## [19] matrixStats_0.63.0          Matrix_1.5-3               
## 
## loaded via a namespace (and not attached):
##   [1] utf8_1.2.3                reticulate_1.26          
##   [3] tidyselect_1.2.0          RSQLite_2.2.20           
##   [5] AnnotationDbi_1.56.2      htmlwidgets_1.6.1        
##   [7] grid_4.1.2                BiocParallel_1.28.3      
##   [9] Rtsne_0.16                munsell_0.5.0            
##  [11] ScaledMatrix_1.2.0        codetools_0.2-18         
##  [13] ica_1.0-3                 statmod_1.5.0            
##  [15] xgboost_1.7.3.1           future_1.31.0            
##  [17] miniUI_0.1.1.1            withr_2.5.0              
##  [19] spatstat.random_3.1-3     colorspace_2.1-0         
##  [21] progressr_0.13.0          filelock_1.0.2           
##  [23] highr_0.10                knitr_1.42               
##  [25] rstudioapi_0.14           ROCR_1.0-11              
##  [27] ggsignif_0.6.4            tensor_1.5               
##  [29] listenv_0.9.0             labeling_0.4.2           
##  [31] GenomeInfoDbData_1.2.7    polyclip_1.10-4          
##  [33] farver_2.1.1              bit64_4.0.5              
##  [35] parallelly_1.34.0         vctrs_0.6.1              
##  [37] generics_0.1.3            xfun_0.37                
##  [39] BiocFileCache_2.2.1       R6_2.5.1                 
##  [41] ggbeeswarm_0.7.1          rsvd_1.0.5               
##  [43] locfit_1.5-9.7            bitops_1.0-7             
##  [45] spatstat.utils_3.0-1      cachem_1.0.6             
##  [47] DelayedArray_0.20.0       assertthat_0.2.1         
##  [49] promises_1.2.0.1          scales_1.2.1             
##  [51] beeswarm_0.4.0            gtable_0.3.3             
##  [53] beachmat_2.10.0           globals_0.16.2           
##  [55] goftest_1.2-3             rlang_1.1.0              
##  [57] splines_4.1.2             rstatix_0.7.2            
##  [59] lazyeval_0.2.2            broom_1.0.4              
##  [61] spatstat.geom_3.0-6       yaml_2.3.7               
##  [63] reshape2_1.4.4            abind_1.4-5              
##  [65] backports_1.4.1           httpuv_1.6.8             
##  [67] tools_4.1.2               ellipsis_0.3.2           
##  [69] spatstat.core_2.4-4       jquerylib_0.1.4          
##  [71] RColorBrewer_1.1-3        ggridges_0.5.4           
##  [73] Rcpp_1.0.10               plyr_1.8.8               
##  [75] sparseMatrixStats_1.6.0   progress_1.2.2           
##  [77] zlibbioc_1.40.0           purrr_1.0.1              
##  [79] RCurl_1.98-1.10           prettyunits_1.1.1        
##  [81] rpart_4.1.19              deldir_1.0-6             
##  [83] pbapply_1.7-0             viridis_0.6.2            
##  [85] cowplot_1.1.1             zoo_1.8-11               
##  [87] ggrepel_0.9.3             cluster_2.1.4            
##  [89] magrittr_2.0.3            data.table_1.14.8        
##  [91] scattermore_0.8           lmtest_0.9-40            
##  [93] RANN_2.6.1                fitdistrplus_1.1-8       
##  [95] hms_1.1.2                 patchwork_1.1.2          
##  [97] mime_0.12                 evaluate_0.20            
##  [99] xtable_1.8-4              XML_3.99-0.13            
## [101] gridExtra_2.3             compiler_4.1.2           
## [103] tibble_3.2.1              KernSmooth_2.23-20       
## [105] crayon_1.5.2              htmltools_0.5.4          
## [107] mgcv_1.8-41               later_1.3.0              
## [109] tidyr_1.3.0               DBI_1.1.3                
## [111] dbplyr_2.3.0              MASS_7.3-58.2            
## [113] rappdirs_0.3.3            car_3.1-1                
## [115] cli_3.6.1                 parallel_4.1.2           
## [117] metapod_1.2.0             igraph_1.3.5             
## [119] pkgconfig_2.0.3           sp_1.6-0                 
## [121] plotly_4.10.1             spatstat.sparse_3.0-0    
## [123] xml2_1.3.3                vipor_0.4.5              
## [125] bslib_0.4.2               dqrng_0.3.0              
## [127] XVector_0.34.0            stringr_1.5.0            
## [129] digest_0.6.31             sctransform_0.3.5        
## [131] RcppAnnoy_0.0.20          spatstat.data_3.0-0      
## [133] Biostrings_2.62.0         rmarkdown_2.20           
## [135] leiden_0.4.3              uwot_0.1.14              
## [137] edgeR_3.36.0              DelayedMatrixStats_1.16.0
## [139] curl_5.0.0                shiny_1.7.4              
## [141] nlme_3.1-162              lifecycle_1.0.3          
## [143] jsonlite_1.8.4            carData_3.0-5            
## [145] BiocNeighbors_1.12.0      viridisLite_0.4.1        
## [147] limma_3.50.3              fansi_1.0.4              
## [149] pillar_1.9.0              lattice_0.20-45          
## [151] KEGGREST_1.34.0           fastmap_1.1.0            
## [153] httr_1.4.4                survival_3.5-0           
## [155] glue_1.6.2                png_0.1-8                
## [157] bluster_1.4.0             bit_4.0.5                
## [159] stringi_1.7.12            sass_0.4.5               
## [161] blob_1.2.3                BiocSingular_1.10.0      
## [163] memoise_2.0.1             dplyr_1.1.1              
## [165] irlba_2.3.5.1             future.apply_1.10.0